
#include <stdio.h>
#include <stdlib.h>
#include <cutil.h>
#include <cuTimer.cu>
#include "gpusort_kernel.cu"

void randomInit(float*,int, int, int);
void printResults(float *, int , int, int);
static void sortMatrix(float *, int, int, int);

//Sorting algorithm for the 6.338J Parallel computing class


int main(int argc, char** argv)
{
	CUT_DEVICE_INIT(argc, argv);
	int dimx= 100;
	int dimy= 16384*2;
	int dimz= 96;


	unsigned int size = dimx * dimy* dimz;
	unsigned int mem_size = sizeof(float) * size;
	float* h_matrix = (float*) malloc(mem_size);

	randomInit(h_matrix, dimx, dimy, dimz);
	sortMatrix(h_matrix, dimx, dimy, dimz);


	CUT_EXIT(argc, argv);
}

void printResults(float* h_matrix, int dimx, int dimy, int dimz)
{
	//PRINT FINAL MATRIX
	bool passed = true;
	for(int k = 0; k < dimz; k++)
	{
		//printf( "\nz= %i\n", k);
		for(int i = 0; i < dimx; i++)
		{
			//printf( "x= %i\n", i);
			int sum=0;
			for(int j = 0; j < dimy; j++)
			{
				//printf( "y= %i ", j);
				//printf( "Value= %f\n", h_matrix[k*dimx*dimy + j*dimx +i]);

				if(j>0)
				{
					if(h_matrix[k*dimx*dimy + j*dimx +i] < h_matrix[k*dimx*dimy + (j-1)*dimx +i])
					{
						passed=false;
					}
				}

				sum=sum +h_matrix[k*dimx*dimy + j*dimx +i];
			}
			//printf( "\nsum= %i ", sum);
		}
	}
	printf( "Test %s\n", passed ? "PASSED" : "FAILED");
}



// Allocates a matrix with random int entries
void randomInit(float * h_matrix, int dimx, int dimy, int dimz)
{
    for (int i = 0; i < dimx*dimy*dimz; ++i)
    {
       h_matrix[i]=rand()/(((unsigned)RAND_MAX + 1) / dimy);
    }

    //Print initial matrix
	for(int k = 0; k <dimz; k++)
	{
		//printf( "z= %i\n ", k);
		for(int i = 0; i < dimx; i++)
		{
			//printf( "x= %i\n ", i);
			int sum=0;
			for(int j = 0; j < dimy; j++)
			{
				//printf( "y= %i ", j);
				//printf( "Value= %f\n", h_matrix[k*dimx*dimy + j*dimx +i]);
				sum=sum +h_matrix[k*dimx*dimy + j*dimx +i];
			}
			//printf( "\nsum= %i ", sum);
		}
	}
}

static void sortMatrix(float * h_matrix, int dimx, int dimy, int dimz)
{

	// allocate host memory for matrices
	unsigned int size = dimx * dimy* dimz;
	unsigned int mem_size = sizeof(float) * size;

	float* ht_matrix = (float*) malloc(mem_size);

	for (int i = 0; i < size; ++i)
	{
	   ht_matrix[i]=-1;
	}


	// COPY THE MATRIX FROM THE HOST TO THE DEVICE
	float* d_matrix;
	float* dt_matrix;
	cudaMalloc((void**) &d_matrix, mem_size);
	cudaMemcpy(d_matrix, h_matrix, mem_size,cudaMemcpyHostToDevice);
	cudaMalloc((void**) &dt_matrix, mem_size);
	cudaMemcpy(dt_matrix, ht_matrix, mem_size,cudaMemcpyHostToDevice);

	cudaError_t error0= cudaGetLastError();
	printf("Errors allocating memory: , %s,\n", cudaGetErrorString(error0));
	printf("Sorting...\n");

	// INVOKE SORTING ALGORITHM ON THE GPU
	int blockY= dimz;
	int blockX= (dimy/ TPB) * dimx;
	dim3 dimBlock(TPB);
	dim3 dimGrid(blockX, blockY);

	float totalTime=0.0;

	cuResetTimer();
	mergesort<<<dimGrid,dimBlock >>> (d_matrix,dimx, dimy,dimz);
	float tA=cuGetTimer();
	cudaThreadSynchronize();
	float tA2=cuGetTimer();
	printf("Cuda kernel Launch =%.1f usec, finished=%.1f msec\n",tA*1.e3,tA2);
	totalTime= totalTime + tA2;

	int strip= TPB;
	int k=2;
	while(strip != dimy)
	{
		strip = TPB *k;
		cuResetTimer();
		float tB=cuGetTimer();
		mergeblocks<<<dimGrid, dimBlock >>> (d_matrix, dt_matrix,k,dimx, dimy,dimz);
		cudaThreadSynchronize();
		swapvectors<<<dimGrid, dimBlock >>> (d_matrix, dt_matrix,dimx, dimy,dimz);
		cudaThreadSynchronize();
		cloneholes<<<dimGrid, dimBlock >>> (d_matrix, dt_matrix,k,dimx, dimy,dimz);
		cudaThreadSynchronize();
		float tB2=cuGetTimer();
		printf("Cuda kernel Launch =%.1f usec, finished=%.1f msec\n",tB*1.e3,tB2);
		totalTime= totalTime + tB2;

		cudaError_t error2= cudaGetLastError();
		printf("Errors in kernels:, %s,\n", cudaGetErrorString(error2));

		k= k*2;

	}

	// COPY THE MATRIX BACK FROM THE DEVICE TO THE HOST
	//CUT_CHECK_ERROR("Kernel execution failed");
	cudaMemcpy(h_matrix, d_matrix, mem_size, cudaMemcpyDeviceToHost);
	cudaFree(d_matrix);
	cudaFree(dt_matrix);

	printResults(h_matrix, dimx, dimy, dimz);

	printf("\nSorting Time =%.1f msec\n",totalTime);


}
